'''
Created on Feb 2, 2010

@author: oabalbin
'''

#import sys;
#import string;
#import StringIO;
import os;
import sys;
import glob;
import time
import MySQLdb;
import numpy as np
from optparse import OptionParser
from datetime import datetime
from datetime import date

# My libraries
import signatures.parsers.read_microarray_data as pm
import signatures.parsers.read_gene_lists as gpm
import signatures.db.tables as dt
import signatures.db.query as dq
import signatures.common.classes as cm #record, arrayMatrix
import signatures.preprocess.tools as tl #filter_constant_genes, get_variable_genes, 
import signatures.seedgenes.seedlist as sd #get_seed_genes, write_varin_file
import signatures.bfrm.paramfile as bp
import signatures.bfrm.bfrm_run as br

class pipeline():
    
    def __init__(self, dbhost, dbuser, dbpasswd, database):
        
        self.mt = dt.tables(dbhost, dbuser, dbpasswd, database)
        self.mq = dq.query(dbhost, dbuser, dbpasswd, database)
        #table information
        self.tableDescription = 'probeID VARCHAR(40), geneSym VARCHAR(40), sample VARCHAR(40), log_fold FLOAT(20,14), INDEX probe (probeID), INDEX gene (geneSym), INDEX sample (sample)'
        self.tableName = 'NCI60'
        
        # Parameters for pre-processing data:
        # percentage of low variable gene to be filter out of the analysis
        self.lvar=25
        # percentage of low median gene to be filter out of the analysis
        self.lmed=25
    
    
    def parse_input_annotfiles(self, annotInputFile):
    
        expr_file, seed_file, outdirs, db_file, samples, bfrm, typefile={},{},{},{},{},{},'nofile'
    
        for line in annotInputFile:
            fields = line.strip('\n')
                 
            if fields[0][0] == '':
                continue
            
            if fields[0][0] == '#':
                typefile = fields
                continue
            
            fields = fields.split('=')
            dict = {'#Array expression':expr_file, '#Seed genes':seed_file,'#Database file':db_file,'#Samples':samples,'#Output Directories':outdirs,\
                    '#BFRM':bfrm}[typefile]
            dict[fields[0]] = fields[1]
            
        return expr_file, seed_file, db_file, samples, bfrm, outdirs
    

    def read_files_folder(self,folderpath,ext):
        ''' '''
        # Read files in folder
        myfiles=[]
        for infile in glob.glob( os.path.join(folderpath, '*'+ext) ):
            myfiles.append(infile)
            #print "current file is: " + infile
        return myfiles
    
    
    def process_microarray_files(self,exp_files,tp):
        dbf=[]
        for f in exp_files:
            of=f+'_dbfile'
            tp.microarray_nci60_file(open(f),open(of,'w'))
            dbf.append(of)
        
        return dbf
    
    def create_database(self,dbfiles,doit=False):
        """
        It creates a databases if doit = True
        """
        if doit:
            print 'Making tables'            
            for dbf in dbfiles:     
                if not(self.mt.table_existence(self.tableName)):
                    self.mt.create_table(self.tableName,self.tableDescription)
                    
                self.mt.load_data_table(self.tableName,dbf)
                
    def check_create_dir( self, full_path_name ):
        """
        creates an output folder time stamped
        """
        t = datetime.now()
        outfolder = full_path_name + t.strftime("%Y_%m_%d_%H_%M")
        if not os.path.isdir(outfolder):
            os.mkdir(outfolder)
        
        return outfolder
     
    def write_pathway_activity(self,pathactivity,samplelist, outfolder):
        
        avgf = open(outfolder+'avg_sign.txt','w')
        medf = open(outfolder+'med_sign.txt','w')
        stdf = open(outfolder+'std_sign.txt','w')
        
        avgf.write('Factor'+'\t'+",".join(samplelist).replace(",", "\t") + '\n')
        medf.write('Factor'+'\t'+",".join(samplelist).replace(",", "\t") + '\n')
        stdf.write('Factor'+'\t'+",".join(samplelist).replace(",", "\t") + '\n')
        
        
        for i, fact in pathactivity.iteritems():
            avgf.write(str(i+1)+'\t'+",".join(map(str,fact[0])).replace(",", "\t") + '\n')
            medf.write(str(i+1)+'\t'+",".join(map(str,fact[1])).replace(",", "\t") + '\n')
            stdf.write(str(i+1)+'\t'+",".join(map(str,fact[2])).replace(",", "\t") + '\n')

##########

if __name__ == '__main__':
    
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-f", "--annotFile", dest="annotFile",
                            help="annotation file for all files to use")
    optionparser.add_option("-d", action="store_true", dest="dbtrue",
                            help="create a database. Bool")

    
    (options, args) = optionparser.parse_args()
       
    # start the pipeline. 
    # open the database
    pipe = pipeline("localhost", "oabalbin", "oscar", "signatures")
    expr_folder, seed_file, db_file, samples, bfrm, outdirs = pipe.parse_input_annotfiles(open(options.annotFile))    
    #read folder for expression profiles
    exp_files = pipe.read_files_folder(expr_folder['expdir'],'.txt')
    
    ####### output files
    
    #outfolder = pipe.check_create_dir( outdirs['outdir']+'bfrm/' )
    
    outfile1 = outdirs['outdir']+'bfrm/dataset_mat.txt'  # Expression values
    outfile2 = outdirs['outdir']+'bfrm/dataset_gen.txt'  # geneNames
    outfile3 = outdirs['outdir']+'bfrm/varin.txt'        # seed genes
    outfile4 = outdirs['outdir']+'bfrm/parameters.txt'   # bfrm parameters file
    
    
    ######### Get raw data  ############
    # Get the database files
    tp = pm.parser()
    dbfiles = pipe.process_microarray_files(exp_files,tp)
    pipe.mt.create_connector()
    
    gtp = gpm.geneparser()
    ######## Create database
    
    pipe.create_database(dbfiles,options.dbtrue)
    
    ########## Make queries ############
    
    pipe.mq.create_connector()
    # list of samples:
    samplelist = tp.list_of_names_in_line(open(samples['samples']))
    # bfrm_factors = list of of factors and weights
    # geneinfactor = list of all genes in each factor
    bfrm_factors, geneinfactor = gtp.list_of_names_in_factor(open(seed_file['allgene']))
    print 'The length of gene in factors is ' + str(len(geneinfactor))
        
    # Do query to obtain the expression matrix related with the gene factors
    # Note: In the query function 0s are mask, in RMA 0 is log2(1), so they could be allowed.
    # The NCI60 exp vals have min 1. So they are not affected. But be aware of the masking
     
    expMat = pipe.mq.get_exp_mat_2(samplelist, geneinfactor, pipe.tableName)
    expMat.create_gen_index()
    
    # Even though NCI60 data is RMA normalized.  
    # Oncomine Normalization is applied. 
    expMat.center_scale_norm()
    
    # Calculate the pathway activity as the weighted average of the expression 
    # value of the genes in the factor.
    pathactivity={}
    for i,fact in enumerate(bfrm_factors):
        wavg, wmed, wstd = expMat.weigthed_exp_aver(fact.genedic)
        pathactivity[i] = [wavg, wmed, wstd]
    
    pipe.write_pathway_activity(pathactivity,samplelist, outdirs['outdir'])

    sys.exit(0)
    ######### Pre-process Raw Data ###############
    # NOTE: This parameter pipe.lvar=25 can be modified to control variability filter
    # Get the index of genes with var >=lvar, and med >=lmed 
    
    #genelist = tl.filter_out_constant_genes(expMat, pipe.lvar, pipe.lmed, True)       
    # Get the index and name of the variable genes for a particular experiment
    #mgenelist, dictmgl =  tl.get_variable_genes(expMat,genelist)
 
    # list of external genes to be used as the "universe of genes"
    # e.i: this could be the ERG interactome
    extgenlist = tp.list_of_names_in_line(open(seed_file['univfile']))
    
    # Index of genes belonging to the "universe" which also have var .=lvar and med >=lmedS
    finalgenes, genenotfound = sd.get_gen_index(expMat,extgenlist)
    
    #print finalgenes
    print str(len(finalgenes)), str(len(genenotfound))
    
    ######## Generate the filtered expression matrix for bfrm
    
    expMat_bfrm = tl.get_filtered_expMat(expMat,finalgenes,len(samplelist))
    
    #print expMat_bfrm.geneSynd
    print str(len(finalgenes)), str(len(genenotfound))
        
    ########## Get seed genes
    seedgenes = tp.list_of_names(open(seed_file['seedfile']))
    #seedgenes = tp.list_of_names_in_line(open(seed_file['seedfile']))
    seeding, seednotfound = sd.get_gen_index(expMat_bfrm,seedgenes)
    #seeding, seednotfound = sd.get_seed_genes(finalgenes,seedgenes)
    print len(seedgenes)
    ######## Write expression matrix for variable genes
    ######## write in varin.txt file the index of seed of genes. To use in BRFM
    
    
    # Choose as seed all genes in the universe
    #sd.write_varin_file(expMat_bfrm.genInd.values(),open(outfile3,'w'))
    #pipe.mq.write_variable_genes(expMat_bfrm, expMat_bfrm.genInd.values(), open(outfile1,'w'), open(outfile2,'w'))
    
    # Choose as seed the genes given by the seedfile
    sd.write_varin_file(seeding.values(),open(outfile3,'w'))
    pipe.mq.write_variable_genes_nfact(expMat_bfrm, seeding.values(), open(outfile1,'w'), open(outfile2,'w'), True)

    #sys.exit(0)
    ######### RUN BFRM #######
    ######### Create the parameters file
    # Call to parameters function
    '''
    parmeter_val(self, NLatentFactors, NObservations, NVariables, DataFile, EvolVarIn, EvolVarInFile,
                  EvolIncludeVariableThreshold=0.75, EvolIncludeFactorThreshold=0.75,
                  EvolMaximumFactors=10, EvolMaximumVariables=150):
    '''
    NLatentFactors = 11
    NObservations = len(samplelist)
    NVariables = len(finalgenes)
    DataFile = outfile1
    EvolVarIn = len(seeding)
    EvolVarInFile = outfile3
    EvolIncludeVariableThreshold=0.75
    EvolIncludeFactorThreshold=0.75
    EvolMaximumFactors=20
    EvolMaximumVariables=200
    
    print EvolVarIn
    print EvolMaximumVariables

    
    bfbp = bp.paramfile()
    
    bfbp.parmeter_val(NLatentFactors, NObservations, NVariables, DataFile, EvolVarIn, EvolVarInFile, 
                      EvolIncludeVariableThreshold, EvolIncludeFactorThreshold, EvolMaximumFactors, EvolMaximumVariables)
    
    bfbp.write_bfrmparam_file(open(outfile4,'w'))
    
    #sys.exit(0)
    ########## Create an bfrm_run object 
    bfrm = br.bfrm_run(bfrm['bfrm_p'],bfrm['bfrm'])
    bfrm.execute_bfrm()
    sys.exit(0)